Anomalous probability of large amplitudes in wave turbulence 
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Time evolution equation for the Probability Distribution Function (PDF) is derived for system 
of weakly interacting waves. It is shown that a steady state for such system may correspond to 
strong intermittency. 



Introduction — Wave Turbulence (WT) is a com- 
mon name for the fields of dispersive waves which are 
engaged in stochastic weakly nonlinear interactions over 
a wide range of scales. Numerous examples of WT are 
found in oceans, atmospheres, plasmas and Bose-Einstein 
condensates [1-9]. For a long time, describing and pre- 
dicting the energy spectra was the only concern in WT 
theory. More recently, some attention was given to the 
study of turbulence intermittency. WT intermittency, 
or "burstiness" of the turbulent signal, was observed 
experimentally and numerically and was attributed, as 
in most turbulent systems, to the presence of coher- 
ent structures. Examples include collapsing filaments 
in Bose-Einstein condensates with attractive potentials 
[9,10], condensate quasi-solitons in systems with repul- 
sive potentials [9,11,12], white caps of sea waves at small 
scales [13], freak ocean waves at larger scales [14]. Often, 
such coherent structures are intense but quite sparse so 
that in most of the space waves remain weakly nonlinear 
and mostly unaffected by these structures. 

Recent analysis of the higher order cumulants [15] 
showed that WT becomes strongly non-Gaussian at the 
same length scale where it fails to be weakly nonlinear. In 
scale invariant systems, the ratio of nonlinear time to the 
linear wave period grows as a power-law either in to small 
or toward large wavenumbers. When this growth coin- 
cides with the cascade direction then one expects the WT 
breakdown if the inertial range is large enough. Other- 
wise intermittency never occurs provided that turbulence 
is weak at the forcing scale [16]. Further, even if a signif- 
icant non-Gaussianity occurs, it does not in itself imply 
intermittency because PDF may remain, in principle, of 
the same order as Gaussian in all of its parts. This mo- 
tivates us to study PDFs in WT. The first study of PDF 
in the WT context was done in [6] for the 3-wave sys- 
tems, whereas here we will be concerned with the 4- wave 
case. We are also motivated by a puzzling numerical ev- 
idence of a low-wavenumber intermittency in the system 
of water-surface gravity waves [17] whereas the analysis 
of [15] predicts intermittency at high wavenumbers only. 
Explaining this fact could shed light on the phenomenon 
of freak waves [14]. 

The idea of the present letter is based on the obser- 
vation that even if the "hard" breakdown (as in [15]) 
does not occur, there always be a part of the PDF tail 
for which the amplitudes are too high for WT to work. 



Such a "mild" breakdown will modify the PDF tail in 
a way that may correspond to intermittency. In fact, 
this case is easier to study analytically because WT still 
works for most of the PDF and the wave breaking phe- 
nomenon can be modeled simply as a phcnomcnological 
cutoff of the PDF tail reflecting the fact that no waves 
exist above the breaking amplitude. The wave breaking 
causes "leakage" and, therefore, a flux in the amplitude 
space which is the key phenomenon leading to deviations 
from the Gaussian equilibrium and intermittency. Note 
an analogy with the well-known k-space fluxes (cascades) 
corresponding to Kolmogorov turbulence which is qual- 
itatively different from the thermodynamic equilibrium 
state. In this paper we will derive an equation for the 
wave amplitude PDF and we will find its steady state 
solutions corresponding to the finite flux in the ampli- 
tude space. Consequently, we will show that the result- 
ing wave fields are intermittent at each wavenumber with 
an anomalously large probability of the large-amplitude 
waves. 

Definition of RPA fields — Previously, the ran- 
dom phase approximation (RPA) has typically assumed 
that the phases evolve much more rapidly than the ampli- 
tudes and, therefore, there exist time intervals where the 
phases are random but the amplitudes are determinis- 
tic [1]. However, numerical simulations indicate that the 
phase and the amplitude vary at the same time scale [19]. 
Thus, we need to generalize RPA to the case where both 
the phases and the amplitudes are random quantities. 
Such generalization was done in [18] where higher mo- 
ments for 3-wave systems were considered. In the present 
letter, we will be dealing with 4- wave systems and we will 
work directly with PDF's rather than moments. 

Let us consider a wavefield a(x, t) in a periodic box of 
volume V and let the Fourier transform of this field be 
afe where k £ Z d and d is the space dimension. Later we 
take the large box limit in order to consider homogeneous 
wave turbulence. Let us write complex ak as <Xk = Akipk 
where Ak € 1Z + is the amplitude and ipk € 5 1 is a phase 
factor (5 1 being the unit circle in the complex plane). 
We say the wavefield ak is of the RPA type if all variables 
in the set {Ak,ipk',k € Z d } are statistically independent 
random variables and i/'fe's are uniformly distributed on 
S 1 . Defined this way RPA refers not only to the phase 
but also the amplitude statistics and therefore we suggest 
a slightly different reading of this acronym: "Random 



Phase and Amplitude" . 

The above properties are sufficient for our WT anal- 
ysis and yet such fields may be strongly non-Gaussian. 
Indeed, RPA allows any shape of the PDF for amplitudes 
Ak and, therefore, it will be a good tool for describing 
intermittency. 

Weakly nonlinear evolution — Consider a weakly 
nonlinear wavcficld dominated by the 4-wave interac- 
tions, e.g. the water-surface gravity waves [1,5,7,13], 
Langmuir waves in plasmas [1,3] and the waves described 
by the nonlinear Schrocdinger equation [9]. In the finite 
box, we have the following Hamiltonian equations for the 
Fourier modes of this field, 

ih=eY J Wllb a b,Ke^>8% (1) 

where bi is the wave action variable in the interaction 
representation, I E Z d , <~ 1 is an interaction coef- 
ficient, (jj 1 ^, = u>i + (J a — ui^ — oj v , u>i is the frequency of 
mode I and 6 < 1 is a nonlinearity parameter. We are 
going expand in e and consider the long-time behavior 
of a wave field, but in order to make such an analysis 
consistent we have to renormalize the frequency of (1) as 



ik = e]T Wl a „a aafl a„e iQ >6 l « - n m , (2) 

ot\xv 

where a x = fo ( e 4f2it , Qi = 2e ^ W^\A^(0)\ 2 is the non- 
linear frequency shift arising from self-interactions and 

w{S - w{S + + n a - % - si v . 

For small nonlinearity, the linear time-scale 2ir/ui is 
a lot less than the nonlinear evolution time which (as 
will be evident below, see e.g. (12)) is 2-Kt 2 /u. Thus, 
to filter out fast oscillations at the wave period, let us 
seek for the solution at an intermediate time T such that 
2tt/lu <C T <C 1/uje 2 . Now let us use a perturbation 
expansion in small e, ai(T) = a|°' + eaj 1 ' + e 2 a!f\ Sub- 
stituting this in (2) we get in the zeroth order a time 
independent result, af\T) = &;(0). For simplicity, we 
will write aj(0) = a;, understanding that a quantity is 
taken at T = if its time argument is not mentioned 
explicitly. The first iteration of (2) gives 

(T) = -i W£a aaii a v 6 l °A% + i% t T. (3) 

where A{£ = A^(T) = (e iQ '"» T - l)/iu l £ v . Iterating one 
more time we get 



+ - e J2 ^iW]Z5%a a a^a v E{^ v ,Q)-W]Z6 1 ^ J re iQ >dr \ , with E(x, y) = A(x — y)e iyt dt. (4) 



Evolution of statistics — We will now develop 
a statistical description via averaging over the initial 
fields cifc(O) which are taken to be of the RPA type. 
Of course, to have a non-trivial description valid over 
the nonlinear evolution time, the fields must remain of 
the RPA type over the nonlinear time in the leading or- 
der in e. The proof of this involves considering the full 
multi-particle PDF and will be published separately be- 
cause it is rather lengthy and outside of the scope of 
the present paper [20]. Let us introduce a generating 
function Z(X,t) = ( e A l a fe(*)l ) ; where A is a real param- 
eter. Then PDF of the wave intensities s — \a,k(t)\ 2 at 
each k can be written as an inverse Laplace transform, 
P(s,t) = (S(\a k (t)\ 2 - S )) =^-J+£z(\,t)e-^d\. For 
the one-point moments we have 

MM = (\a k \ 2 n = (\a\ 2p e^ 2 )\ x = = 

rOO 

Z x ... x \\=o = / S p P(s,t)d s , (5) 
Jo 

where p € M and subscript A means differentiation with 



respect to A p times. 
At t = T we have 

Z(T) = (e^+^M 2 '! 2 ) = 

( e ^t 0, l 2 (l + A e(a «4 0) +cc) 
Xe 2 (\a^\ 2 + (a^ + cc)) + 

^f(a^ +cc)%) A 

= Z(0) + eA( e A l a - 0) l 2 (4 1) 4 0) +cc)^) A + 
e 2 (((X + X 2 A 2 )\a^\ 2 + X(a^a^+cc) 

+ y(4 1)2 4° )2 + cc)M.4 (6) 

where cc stands for complex conjugate of the previous 
terms and (. . .),/, and (. . .) a denote phase and ampli- 
tude averaging respectively. Note that in RPA fields the 
phases and the amplitudes are statistically independent 
so that these two averaging could be done independently. 



First let us substitute and aj^ from (3) and (4) 
respectively and perform the phase averaging. For the 
terms proportional to e we have 



-2iJ2wtA 2 k Al.T + i^AlT. (7) 



where we have used the fact that A(0) = T and we have 
used the RPA's "Wick's Theorem" 

(a k a a a> l a„U=AlA 2 a (6ffi + 6 k JZ). 

Note that the above expression should also contain the 
singular cumulant, i.e. term (A* — 2A 2 t )S"8^S k , see [18]. 
However we do not write this term here, since its contri- 
bution is of the order of 1/N smaller because it has one 
extra delta function. 

We see from (7) that the choice 



n k = 2^2w^Al 



(8) 



makes the contribution of (4^ terms to be equal 

to zero. 

We therefore obtain 

Z{T) - Z(0) = 

A £ 2 ((Ai4 l )| 2 (i + A| flfc | 2 ) + 4 2 )4°) + 4 2) 4 0) + 

+ ^((4° ) 4 1) ) 2 + cc)),) A (9) 
where as an intermediate result we have 

<4 0) 4 2) >. = 

2 £ SfrlWfr \ 2 A 2 k {A\Al 2A 2 a Al)E(0, Co k «) 



and 



a\\2 



Here terms proportional to T 2 drop from (4°^ 4^)^ an( l 
(|4^I 2 ) because of the choice (8) of frequency renormal- 
ization. Furthermore, 

((4 1) 4 0) ) 2 >, = 

-2W kk A e k T 2 n k - A\Vl\T 2 -^ W ^A k A a Tf 

a 

To complete the derivation of the equation for the time 
evolution of the generating function Z(T) we have to 
take a large box limit, which implies that sums will be 



replaced with integrals, the Kronecker deltas will be re- 
placed with Dirac's deltas, 5 l £ n — ► S!^ n /V, where we in- 
troduced short-hand notation, 5%~ n — S(k a +ki — k m — k n ). 
Then (9) will still hold, but with 

(4 0) 4 2) ). = 

2 J dl23S^\W k i\ 2 Al(A 2 A 2 -2A 2 A 2 )E(0^ kl 3 ) 



and 



|4 1) | 2 ) V , = J di236%\w£\ 2 A 1 A 2 A 3 \A(u>%)\ 2 . 



We also have 



{( a k a k ) /</> - u > 



(10) 



because this terms will be 1/N times smaller than 
(|a^| 2 )^, and (a k a k 2 ^}^, terms because it has one less 
summation index. Therefore it vanishes in the N — > oo 
limit. 

Further we take a large T limit, and take into account 
that 

.1. 



and 



lim E(0,x) = 0k6(x) + iP(-))T, 

T-»oo X 



lim |A(a;)| 2 = 2nT6(x), 

T — >oo 



(see e.g. [2]). 

Finally we perform amplitude averaging, noticing that 



Z(0) = (e A|afcl ) A , 



and 



8Z_ 

ax 



= \\a k \ e 



2 A|aJ 2 \ 

1 ]A- 



to obtain 

Z{T) = Z(0) + e 2 T ■ (X V Z + (\ 2 r, - A 7 )Z A ). 
Approximating (Z(T) - Z(0))/T by Z, we have 
Z = XvZ + (X 2 T]- Xn-f)Z x , 

where 

r} k =A-ne 2 J \W^\ 2 5^6(w^) ni n 2 n3dl23, 

7fe = 8tt£ 2 J \ W^\ 2 5%l5{u k l) m(n 2 + n 3 ) - n 2 n 3 



(11) 



(12) 



dl23, 



here wavenumbers k, k\, k 2 , k 3 € TZ d , 5's now mean Dirac 
5-functions, ni,2,3 = ^(^1,2,3) and dl23 = dkirfk 2 rfk 3 , 



Differentiating (12) with respect to A p times we get the 
evolution equation for the moments: 

Mi p) =-p lk M { k p) +p 2 Vk M[ p - 1 \ 

which, for p = 1 gives the standard kinetic equation, 
rife = — 7ferife + rjh- First-order PDE (12) can be easily 
solved by the method of characteristics. Its steady state 
solution is 

Z=(l-Xn k )' 1 

which corresponds to the Gaussian values of momenta 
MW = p\rP k . However, these solutions are invalid at 
small A and high p's because large amplitudes s = \a\ 2 , 
for which nonlincarity is not weak, strongly contribute in 
these cases. Due to the integral nature of definitions of 
and Z with respect to the s = \a\ 2 , the ranges of 
amplitudes where WT is applicable are mixed with, and 
contaminated by, the regions where WT fails. Thus, to 
clearly separate these regions it is better to work with 
quantities which are local in s = \a\ 2 , in particular the 
probability distribution P(s). Taking the inverse Laplace 
transform of (12) we have 

P + d s F = 0, (13) 

where F = —s(jP + T]d s P) is a probability flux in the 
s-space. Consider the steady state solutions, P = 0, 

-s(fP + r]d s P) = F = const. (14) 

Note that in the steady state 7/77 = n k which follows 
from kinetic equation. The general solution to (14) is 

P Pfiom ~t~ Ppart 

where 

Phom = const exp (s/n) 

is the general solution to the homogeneous equation (cor- 
responding to F = 0) and P par t is a particular solution, 

Ppart = -(F/r])Ei(s/n) exp (-s/n) 

where Ei(x) is the integral exponential function. 

At the tail of the PDF, s ^> n k , the solution can be rep- 
resented as series in 1/s, P par t = —F/{s~f) — r)F/(^s) 2 + 
■ ■ ■ . Thus, the leading order asymptotic of the finite-flux 
solution is 1/s which describes strong intcrmittency. 

Note that if the weakly nonlinearity assumption was 
valid uniformly to s = oo then we had to put F = to en- 
sure positivity of P and the convergence of its normaliza- 
tion, J Pds = 1. In this case P — Phom = n exp (—s/n) 



which is a pure Rayleigh distribution corresponding to 
the Gaussian wave field. However, WT approach fails 
for the amplitudes s > s n i for which the nonlinear time 
is of the same order or less than the linear wave period 
and, therefore, we can expect a cut-off of P(s) at s = s n i. 
An estimate based on the dynamical equation (1) gives 1 
s n i = oj/eWk 2 . This phcnomcnological cutoff can be 
viewed as a wave breaking process which does not allow 
wave amplitudes to exceed their critical value, P(s) = 
for s > s n i- Now the normalization condition can be 
satisfied for the finite-flux solutions. However, having a 
constant negative flux F < corresponds to a source at 
s = s n i which dictates the necessity of a sink for some 
s < s n i to preserve the normalization of P(s). Note how- 
ever that the probability sink does not have to correspond 
to any physical "removal" of waves with certain ampli- 
tudes. The sink should be present solely because the 
probability is diluted due to acceptance of new members 
with s = s n i into the statistical ensemble. In this case, 
the sink must be proportionate to the probability and, 
taking into account the normalization condition, we can 
write a modified equation for the PDF in the presence of 
cutoff, 

P - a s (s 7 p + s V d s p) = —Ft, (15) 

with F„ = —P{s n i)l I s n i- The general solution so- 
lution to this equation is P = [C — F*Ei(s/n — 
log s)/rj\ exp (—s/n), there constant C is fixed by the 
normalization condition. This solution is close to the 
Rayleigh distribution in the PDF core, s ~ n, and it has 
a 1/s tail atn<s<s„i. 

Discussion — We found that the WT intermittency 
shows as an anomalously high (~ 1/s) probability of the 
large-amplitude waves whereas at lower amplitudes dis- 
tribution appears to be close to Rayleigh (~ e~ s /™) which 
corresponds to Gaussian wave fields. We showed that 
wave breaking is essential for WT intcrmittency to be 
present in the system, yet the details of wave breaking 
are not important. The role of wave breaking is just to 
ensure that no wave can have amplitude greater than 
critical value s n i- This simple condition leads to huge 
mathematical consequences as it generates the flux so- 
lutions in the amplitude space and therefore creates the 
1/s intermittency. On the other hand, the amplitude of 
the 1/s tail is not prescribed by WT and will depend on a 
particular wave breaking mechanisms in a particular sys- 
tem. However, some conclusions about the dependence 



1 This estimate assumes that if the wave amplitude at some k happened to be of the critical value s n i then it will also be of 
similar value for a range of fc's of width k. In other words, strong nonlinearity widens the fc-space correlation from zero (RPA 
value) to k (value for the coherent structures involved in the wave breaking). 



of the tail amplitude on the physical parameters can be 
reached using a dimensional arguments. 

Consider a classical example of the gravity waves on 
surface of deep water. The linear dispersion relation is 
given by uik — Vgk, and the coefficient of nonlinear inter- 
action W*" is given in [1]. This system has two power- 
law steady state solutions. First one is the spectrum 
corresponding to the direct cascade of energy toward 
high-wave numbers, oc k~ A [1,4]. Second one is the 
spectrum corresponding to the inverse cascade of wave 
action toward the small k values, n k oc fc~ 23 / 6 . In addi- 
tion to the gravity constant g, the only quantity which 
determines the state of the system in the direct cascade 
range is the energy flux V whereas in the inverse cas- 
cade range - the particle flux Q. The PDF tail strength 
can be characterized by its area which is a dimensionless 
number and, therefore, has to depend on the relevant 
dimensionless combinations in the direct and the inverse 
cascade ranges, V 2 l 3 k x l 3 j g and Qk/g respectively. Thus, 
the PDF tail thickness grows with k but its length de- 
screases until it completely disappears at k <~ k ni (equal 
to g 3 /V 2 and g/Q respectively). 

This effect is illustrated in Figure 1 which shows PDF's 
obtained by numerical simulations of the surface waves 
on deep water forced at low k's and dissipated at high 
fc's. Pseudospectral numerical method similar to that of 
[17], [7] was used on a 256x256 grid. 

At moderate wavenumber (k = 15k min ) one can see 
a PDF tail in the range 4n k < s < lOn^ characterized 
by an order of magnitude enhanced probabilities with 
respect to the Rayleigh distribution. Unfortunately the 
range of s where PDF converged to a stable value in this 
experiment was not large enough to reach s 3> n values 
and, therefore, for an asymptotic scaling to develop. To 
increase this range a much longer computing to gain good 
statistics of very rare events at the PDF tail is necessary, 
which we can not perform with our resources. 

At a higher wavenumber {k = 35fc m i„) one can see that 
the large amplitude waves are less probable than the ones 
predicted by the Rayleigh distribution. This is because 
the wavebreaking happens now closer to the PDF core 
causing the PDF cut-off seen at the figure. 

In this letter we considered WT which is weak on av- 
erage so that the wave breaking occurs only in the PDF 
tail, i.e. s n i 3> n. It does not apply to the cases when, at 
some large fc, the wave breaking may become so strong 
that it occurs for most of the waves in the PDF core. 
These cases where predicted and discussed in [15], but 
their statistics would be hard to describe analytically be- 
cause of the strong nonlinearity. 
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FIG. 1. PDF of |afe| 2 for k = 15k m i„ (thick curve) and 
k = 35fe m i n (thin curve) and their comparison with Rayleigh dis- 
tribution (dotted line). Amplitude s is normalised so that the two 
curves have the same slope at s = 0. 
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